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ABSTRACT 

Context. The understanding of the formation process of massive stars (^8 M Q ) is limited, due to a combination of theoretical com- 
plications and observational challenges. The high UV luminosities of massive stars give rise to chemical complexity in their natal 
molecular clouds, and affect the dynamical properties of their circumstellar envelopes. 

Aims. We investigate the physical structure of the large-scale (~10 4 -10 5 AU) molecular envelope of the high-mass protostar 
AFGL2591. 

Methods. Observational constraints are provided by spectral imaging in the 330-373 GHz regime from the JCMT Spectral Legacy 
Survey and its high frequency extension. While the majority of the ~160 spectral features from the survey cube are spatially unre- 
solved, this paper uses the 35 that are significantly extended in the spatial directions. For these features we present integrated intensity 
maps and velocity maps. The observed spatial distributions of a selection of six species are compared with radiative transfer models 
based on (i) a static spherically symmetric structure, (ii) a dynamic spherical structure, and (iii) a static flattened structure. 
Results. The maps of CO and its isotopic variations exhibit elongated geometries on scales of ~100", and smaller scale substructure 
is found in maps of N 2 H + , o-H 2 CO, CS, S0 2 , C 2 H, and various CH 3 OH lines. In addition, a line of sight velocity gradient is apparent 
in maps of all molecular lines presented here, except SO, S0 2 , and H 2 CO. We find two emission peaks in warm (E up ~ 200 K) 
CH 3 OH separated by 12" (12 000AU), indicative of a secondary heating source in the envelope. The spherical models are able to 
explain the distribution of emission for the optically thin H 13 CO + and C 34 S, but not for the optically thick HCN, HCO + , and CS, nor 
for the optically thin C 17 0. The introduction of velocity structure mitigates the optical depth effects, but does not fully explain the 
observations, especially in the spectral dimension. A static flattened envelope viewed at a small inclination angle does slightly better. 
Conclusions. Based on radiative transfer modeling, we conclude that a geometry of the envelope other than an isotropic static 
sphere is needed to circumvent line optical depth effects. We propose that this could be achieved in circumstellar envelope models 
with an outflow cavity and/or inhomogeneous structure at scales < 10 4 AU. The picture of inhomogeneity is supported by observed 
substructure in at least six different species. 
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1. Introduction 

A collapsing protostellar cloud is heated by gravitational con- 
traction of the envelope and radiation from the central protostar. 
Therefore, the process of star formation depends crucially on the 
disposal of thermal energy through (molecular) transition lines. 
Molecular signatures from star-forming regions are an impor- 
tant probe in studies of star formation. While the energy budget 
is reasonably well understood for protostars with masses < 8 M© 
(e.g., [Evans 1999; Larson 2003 ), the problem is more complex 
for high-mass stars (> 8 M Q ). Not only do they require more ma- 
terial and involve (orders of magnitude) more energy, high-mass 
stars form in more deeply embedded natal clouds, and gener- 
ally in more clustered environments. In addition, high-mass pro- 
tostars are rarely observed due to the fact that they form less 
often, and that if they do, their formation timesc ales are much 
shorter than for their low-mass counterparts (see Zinnecker & 
Yorke 2007 1 for a review). In addition, young high-mass stars 
have the ability to significantly alter the chemical composition 
of their parental clouds by the copious amounts of UV-radiation 



they emit ( Stau ber et al. 2004 ). It is exactly this impact that high- 
mass stars have on the surrounding interstellar medium - both 
when they form and when they explode as supernovae - which 
makes it important to study their formation process in a broader 
context. 

Due to the above challenges in the study of high-mass star 
formation, it is as yet unclear whether a high-mass star forms fol- 
lowing a process similar to low-mass star formation (core con- 
traction — > envelope collapse — > disk accretion — > disk evapo- 
ration), or that other processes such as competitive accretion or 
stellar mergers are at play (e. g.,|Bonnell et al.|2001| Bonnell & 
Bate 2005} |Hocuk &~S paans 2010). In any case, much higher 
accretion rates and more energetic feedback processes must be 
involved in high-mass star formation, providing additional chal- 
lenges in theoretical modeling. Investigations into the physical 
and chemical structure (at various scales) of the molecular en- 
velopes of massive young stars can shed light on the accretion 
process and the energy distr ibution and dissipation in r egions 
of high-mass star formation (Eva ns 1 1999} Cesaroni|2005 ). One 
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hypothesis to consider is whether high-mass star formation in 
its first stages proceeds less isotropically than low-mass star for- 
mation, undergoing non- spherical accretion even before material 
reaches a potential circumstellar disk. 

This paper studies the molecular envelope of AFGL2591, an 
isolated massive star-forming region located in the Cygnus-X re- 
gion at galactic coordinates (£, b) = (78?9, (f.ll). The distance to 



AFGL2591 is between 0.5 and 2.0 kpc (see |Van der Tak et al. 
1999 [[Schneider et al.|2006| ). While the commonly assume d av 



erage distance for the Cygnus-X complex is 1.7 kpc ( [Motte et aL 
2007), the distance to individual objects in this local spiral arm 



structure varies. For consistency with earlier work on this source, 
we adopt a distance of 1.0 kpc and we refer to discussion in Van 



der Tak et al. ( 1999 ) concerning the effects of a larger assumed 
source distanc e. The protostar is estimated by |Van der Tak & 
Menten (2005) to have a mass of 16 M©. The object emits most 



of its energy in the infrared, due to re-radiation of optical/UV ra- 
diation from the central object by the surrounding gas and dust. 
It has a bolometric luminosity of ~2xl0 4 L© (at 1 kpc) an d ex- 
hibits powerful large (> Y) bipolar molecular outflows (Lada 
|et al.|1984[|Hasegawa & Mitchell|1995| ). 

Previous modeling of the molecular envelope of AFGL2591 
has assumed a spherical morphology and found that observations 
of molecular lines and dust can be explained using a power law 



density profile (Van der Tak et al. 1999, 2000, Doty et al. 



Stau ber et al |2005[ |Benz et al.||2007UDe Wit et al.|2009| ). |Van| 



2002; 



der Tak et al.| (| 1999) suggested the presence of cavities in the en 



velope, carved out by the outflows, as a solution to optical depth 
effects which prevent molecular line radiation from escaping the 
cloud if the envelope is spherical and isotropic. After evidence 
of such cavities was observed at high spatial resolution (~ 2") in 
the near- infrared by|Preibisch et aL] ( |2003| ), the idea was invoked 
again by Brude rer et al. ( 2010a|b| ) to allow far-UV radiation from 
the protostar to directly excite molecular species in the envelope. 
The overall structure of the envelope, however, is still assumed 
to be spherical. 

This paper is structured as follows. In Sect. [2] we describe the 
unbiased sub-millimeter spectral imaging survey of the molec- 
ular envelope of AFGL2591; the resulting maps of integrated 
intensity and velocity structure for 35 molecular transitions are 
presented in Sect. [3] While this paper focuses on the spatial do- 
main in order to investigate the large-scale (> 10 4 AU) envelope, 
a full exploitation of the spectral domain will be addressed in 
Van der Wiel et al. (2011b, in prep.). Sect. [4] describes various 
approaches to modeling a selection of six molecular emission 
lines using a radiative transfer code. Section[5j finally, is devoted 
to the discussion of the results of the radiative transfer modeling, 
a phenomenological interpretation of the lines that are not mod- 
eled, and general conclusions. 

2. Observations 

2.1. HARP-B and the JCMT Spectral Legacy Survey 

The observations presented here were taken with the 16-element 
Heterodyne Array Receiver Programme B (HARP-B) and the 
Auto-Correlation Spectral Imaging System (ACSIS) correlator 
( [Smith et al.] [20081 pent et al.|[j500l |Buckle et aLl[2009] ) at 
the James Clerk Maxwell Telescop^(JCMT) on Mauna Kea, 
Hawai'i. The strength of the multipixel HARP-B instrument lies 



1 The James Clerk Maxwell Telescope is operated by the Joint 
Astronomy Centre on behalf of the Science and Technology Facilities 
Council of the United Kingdom, the Netherlands Organisation for 
Scientific Research, and the National Research Council of Canada. 



in the combination of high-resolution (1 MHz, ^1 kms *) het- 
erodyne spectroscopy with instantaneous mapping capabilities. 

The JCMT Spectral Legacy Survey (SLS, |Plume et al.|2007| ) 
has been allocated a total of 456 hours to map four Galactic 
targets in the spectral region between 330 and 360 GHz. The 
targets are: the Orion Bar, a prototypical photodissociation re- 
gion ( [Van der Wiel et al.|[2U09l ); NGC1333 IRAS 4, a low- 
mass star-forming region; W49A, an active star formation com- 
plex ( Roberts et al.|[20TT] ); and the isolated high-mass star- 
forming region AFGL2591, which is the topic of this study. The 
SLS is complemented at higher frequencies (360-373 GHz) by 
HARP maps from separate proposals between 2007 and 2010. 
Observations for the SLS (330-360 GHz) started immediately 
after the commissioning of the HARP-B instrument; AFGL2591 
was observed in November 2007, March, July and September 
2008, May-July 2009, and May-July 2010. The high frequency 
maps were observed in August and September 2007 and July 
2008. AFGL2591 is the first of the four SLS targets for which the 
spectral imaging has been completed in the 330-360 and 360- 
373 GHz spectral regions. 

The 16 receptors of the HARP receiver are distributed in a 
4x4 pattern, regularly spaced at 30" intervals. The harp4_mc 
jiggle position switch mode was used to create maps of a 2' x 2' 
area with pixels spaced by 7.5", roughly half of the JCMT 
beam width at the relevant frequencies (15" at 345 GHz). While 
the SLS observations were done in relatively wet weather at 
0.12 < T225GHZ < 0.20, the observations above 360 GHz re- 
quire T225GHZ < 0.08, since the Earth's atmosphere becomes in- 
creasingly opaque as the frequency approaches 375 GHz. Apart 
from the weather constraints, the observing strategy and reduc- 
tion method is identical for the high frequency data and the SLS 
data. 

The spectra of AFGL2591 are calibrated by measurements 
at an off position 340" east of the target. After inspection of 
the calibrated data product, we note that there is a possibility 
of 12 CO 3-2 contaminating emission at the off position. There 
is no evidence of contamination from the off position for other 
spectral lines. Redundancy for bad receptors is built in by con- 
secutive observations at an offset of 15", or by rotation of the 
entire array during periods where more than three of the six- 
teen receptors were not functioning. Redundancy for frequency 
spikes is provided by offsetting the central frequency of consec- 
utive observations, each spanning 1 GHz bandwidth, by 200 or 
400 MHz. Pointing was checked regularly and corrections were 
generally < 2-3". 

2.2. HARP-B data processing 

The raw time series files are processed in blocks of 1 GHz band- 
width, using standard procedures from the Star link software 
package. Each scan file is inspected for bad receptors, base- 
line issues and frequency spikes, which are masked before scan 
files at the same frequency are summed and converted from time 
series into three-dimensional data cubes (RAxDecxfrequency). 
The generally noisy band edges are discarded, typically leav- 
ing 0.82 GHz of usable bandwidth per block. Baselines are fit- 
ted to line-free regions in every cube using a spline algorithm; 
the generally nearly linear baselines are subsequently subtracted. 
The baseline corrected cubes are mosaicked into a larger cube 
spanning (2') 2 on the sky and 330.0-373.3 GHz in frequency. 
The noise level is reduced by exploiting the frequency overlap 
between subsequent frequency settings. The spectrum with na- 
tive frequency bins of 1 MHz (0.8-0.9 kms -1 between 373 and 
330 GHz) is resampled to uniform 1 kms -1 velocity resolution. 
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Fig. 1. Root-mean- square (rms) noise statistics of the data prod- 
uct, calculated in 1 km s" 1 channels. The square symbols denote 
the median rms in all map pixels over a selected frequency range, 
with their error bars indicating the rms of the rms noise across 
all map pixels in the same frequency range. The triangle sym- 
bols mark the rms value in the least noisy (pointing up) and most 
noisy pixel (pointing down) in the map. The horizontal dashed 
line indicates the survey target noise for this object. 



Intensities are corrected for the JCMT main beam efficiency of 
0.61+0.03 ( [Buckle et al.|2009| ) in the 345 GHz window and are 
therefore given in terms of main beam temperature T m \>. To con- 
vert observed T A * values to flux density units (Jy), multiplication 
by a factor 29.4 would be necessary. 

Figure [T] shows the median of the root-mean- square (rms) 
noise in the resulting data cube in some selected line-free fre- 
quency ranges, as well as the spread in rms across the map pix- 
els. The noise in T m \> units in both the 330-360 GHz and the 
360-373 GHz range are generally close to the r m b target noise 
level of 65 mK (converted from 25 mK in 2.5 km s" 1 channels in 
T A * units). The lowest noise values are usually found near the 
spatial center of the cube. 

2.3. Submillimeter continuum images 

The analysis in this paper uses subsidiary submillimeter con- 
tinuum images of AFGL2591. These 450 and 850 /mi im- 
ages have previously been presented by |Van der Tak et aL 
( 2000). The raw images, available from the JCMT archive at the 
Canadian Astronomy Data Centre, are processed using the stan- 
dard ORAC-DR pipeline for SCUBA jiggle-maps. An extinction 
correction is applied, based on sky-dip measurements of the sky 
optical depth. The resulting maps have units of mJy beam -1 and 
are sampled in 3" pixels. The FWHM beam size is 8" for the 
450 yum map, and 14" for the 850 jim map. This makes the latter 
map especially suitable for comparison with our HARP-B spec- 
tral maps. 



3. Observed spatial distribution 

3.1. Selection of spatially resolved maps 

While the majority of the -160 spectral lines show a spatial dis- 
tribution that is confined to the beam size of the JCMT (14- 
15" FWHM at the relevant frequency), in this paper we se- 



lect the lines that are spatially resolved. The selection is made 
by fitting a two-dimensional Gaussian profile to a spatially re- 
sampled (3"75 pixels) integrated line intensity ( J r m bdV) map 
for each spectral feature. Features are registered in Table [T] if 
^fwhm > #fwhm + 3cr fl , with <2fwhm an d °~a me fitted long axis 
and associated error, and #fwhm the JCMT beam size. Spatial 
extent parameters in Table [T] are given as measured in the ob- 
served maps, without beam deconvolution. We select a total of 
35 spectral features with this method. 

The association of spectral features to molecular transitions 
is done using CDMS (Mulle r~et al.|2005| ). In cases where mul- 
tiple molecular transitions fall within a few MHz of a measured 
line frequency, we choose the simplest molecule and/or the low- 
est upper level energy. For the 35 lines selected above, which 
are relatively strong ( J r m bdV > 2Kkms _1 ), the identifications 
listed in Table [T] are unambiguous. Some features are marked as 
a blend of several (hyperfine) transitions of the same or different 
molecules. 

Maps of integrated line intensity for each of the lines with 
extended emission are shown in Figs.[2}j4] Intensity is integrated 
over the velocity interval indicated in Table [2] for each map. 
The intervals are determined manually, with the aim of incor- 
porating all emission belonging to the spectral line in question 
without including any signal from neighboring lines. Intervals 
are generally symmetric around the systemic Vi sr of -5.8 km s" 1 
for AFGL2591, the interval width depending on the width of 
the line (CO 3-2 integrated over [-29, 6] km s" 1 is the most ex- 
treme case). In the maps of integrated line intensity (Figs. [2]- 
[4]), contours are drawn based on resampled maps with pixels 
of 3"75. Contour levels are scaled to the maximum intensity in 
each map listed in Table [2] For comparison, SCUBA 850 /mi 
continuum emission is shown on the background of each inte- 
grated line intensity map. In addition, we show an example of 
a map which is not extended according to our criterion: SO2 
152,14-141,13 in Fig. [5] In any further discussion of SO2 mor- 
phology, we do not consider the 152,i4-14i f i3 mentioned here, 
but only the three transitions 5^-^2,2, 233 ? 2i-232,22, and 232,22- 
23i ? 23 that are listed in Table [T] 

3.2. Spatial distribution of integrated line intensity 

The emission of all molecular tracers from Table [TJ with the ex- 
ception of CO, its isotopologs, and SO2, is distributed over a 
central region -20-50" in diameter (Figs. [2] [3] [4]). The flattened 
morphologies of roughly half of the maps are consistent with a 
position angle of the major axis between 65° and 95° (Table [T]). 
The only maps in which more extended emission is visible are 
those of the 3-2 transition of CO, 13 CO and C 17 (Fig.||, while 
SO2 emission is the most concentrated. Maps of CO and 13 CO 
trace not only the smaller scale morphology already noted above 
for other molecules, but also an additional larger scale envelope 
(up to -100") at lower intensity with a different orientation at a 
position angle of -20°. The large scale (0.07-0.30 pc from the 
ce ntral position) extends to th e region where previous modeling 
by Van der Tak et al. (|2000| ) suggests the molecular hydrogen 
density to be <10 2 cm J and dust and gas temperatures <40 K. 

A 'plume' extending to the north at position angle -10° 
is apparent in maps of CO, C 17 0, HCN, HNC, N 2 H + , HCO + , 
o-H 2 CO, and CH3OH li-Oo A + , 7 -6 A + and 7_i-6_i E 
(Figs. [2j |3j [4]). Its intensity level is generally at -10% of the 
maximum, but it is stronger in the CH 3 OH lines relative to the 
peak intensity. The feature appears to emerge from a position 
roughly 15" east of the peak position. Alternatively, it might 
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Fig. 2. Maps of emission from the 3-2 transitions of CO, 13 CO, C 17 0, transitions of N-bearing molecules, HCO + and H 13 CO + 4- 
3. White contours represent integrated line emission and are at 90%, 80%, . . . , 10% of the maximum intensity, listed in Table 2] 
The contour at 50% is thicker. An extra contour is added for HNC 4-3 at 5% of the maximum intensity. The beam size at the 
relevant frequency is indicated in the top right corner of each panel. To avoid noise in the figures, contours below 0.7 K km s" 1 are 
not shown. SCUBA 850 /mi continuum emission is shown in color scale and thin black contours; the color scale is stretched from 
10 to 5600 mJy beam" 1 . 



trace a warped structure of the collapsing molecular cloud. The 
plume is also apparent in the 850 /mi continuum map (Fig. [2]). 
W e recognize the same featur e in the 13 CO 3-2 map in Fig. 3 
of Van der Tak et al.] ( |1999| ). The northern warped plume is 
not present in maps of C2H nor in any of the sulfur bearing 
molecules. This difference in observed morphology is not just 
an effect of signal-to-noise, since some of the latter transitions 
have peak line strengths (Table [2]) comparable to the maps that 
do show the pl ume . The velocity structure of this plume is de- 
stribed in Sect. 13.31 



Maps of CO and isotopologs, CN 3 V2 -2 5/2 AF=l (Fig. [2]), 
CS, o-H 2 CO (Fig. [3]), C 2 H and CH 3 OH 4 , 3 -3_i, 3 (Fig. [4]) show 
a hint of two separate southern 'arms' ~30" south of their 
main peak, separated by > 20" in the east-west direction. 
Alternatively, this feature can be described as one indentation. 

The spatial distribution of N 2 H + 4-3 emission (Fig. [2] top 
right panel) is strikingly different from those of other molecu- 
lar transitions observed here. It is stretched significantly in the 
east-west direction, mostly due to the spatial separation of the 
velocity components described in Sect. |3.3| Whereas various 
other species manifest two southern 'arms' (see above), only 




28s 



24s 20h29m20s 



RAQ2000) 

Fig. 5. Spatial distribution of integrated S0 2 152,i4-14i i 3 emis- 
sion (white contours), which is not extended. The colorscale and 
black contours represent 850 /mi continuum, like in Fig. [2] 
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Fig. 3. Like Fig. [2] but for H2CO, SO, SO2, and CS. Maximum intensities are listed in Table^ 
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Fig. 4. Like Fig. [2J but for CH3OH and C2H. Maximum intensities are listed in Table [2] 
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Table 1. Molecular lines with observed spatially extended emission. 



Molecule 


Transition 


Frequency 61 


ARA b 


ADed 7 


Angular extent (Gaussian FWHM) C 


Position angle^ 






(MHz) 


(") 


n 


Major axis (") 


Minor axis (") 


(degrees) 


CO 


7 = 3-2 


345796.0 


1.0 + 0.2 


4.7 + 0.4 


84.6 + 1.6 


53.9 + 0.6 


94.65 + 0.02 


13 co 


7 = 3-2 


330588.0 


2.7 + 0.2 


-1.4 + 0.2 


53.9 + 0.6 


46.6 + 0.5 


121.88 + 0.05 


c 17 o 


7 = 3-2 


337061.2 e 


-2.9 + 0.2 


-2.6 + 0.2 


42.4 + 0.6 


30.5 + 0.4 


78.05 + 0.03 


C 2 H 


Nj = 4 9/2 - 3 7/2 


349338.3 / 


-2.5 + 0.2 


-1.3 + 0.1 


30.7 + 0.5 


22.0 + 0.3 


105.91 + 0.03 


C 2 H 


Nj = 4 7/2 - 3 5/2 


349340.0* 


-2.8 + 0.2 


-0.9 + 0.1 


29.1 +0.4 


20.6 + 0.3 


108.54 + 0.03 


CH3OH 


J K = 7i - 61 A + 


335582.0 


0.2 + 0.4 


-3.4 + 0.4 


23.8 + 1.1 


13.7 + 0.6 


135.41 + 0.06 


CH3OH 


J K =l2 l -12 A T 


336865.1 


-1.7 + 0.7 


-1.7 + 0.6 


21.2+ 1.8 


13.1 + 1.1 


118.79 + 0.12 


CH3OH 


Jk = 7 - 6 ( , E 


338124.5 


1.3 + 0.3 


-2.6 + 0.3 


21.5 + 0.8 


17.5 + 0.7 


34.25 + 0.13 


CH3OH 


J K = l- l -6. 1 E 


338344.6 


3.2 + 0.4 


-2.2 + 0.3 


30.3 + 1.0 


20.6 + 0.7 


89.08 + 0.06 


CH3OH 


J K = l -6 A + 


338408.7 


4.7 + 0.4 


-2.3 + 0.3 


31.2 + 0.9 


23.0 + 0.7 


70.66 + 0.07 


CH3OH 


^ = 7 2 -6 2 E h 


338722.3^ 


0.7 + 0.3 


-3.2 + 0.3 


19.6 + 0.8 


15.5 + 0.6 


114.61 +0.12 


CH3OH 


Jk = 4 - 3_! £ 


350687.7 


1.8 + 0.3 


-1.4 + 0.3 


34.6 + 0.8 


25.3 + 0.6 


72.59 + 0.05 


CH3OH 


7 i , = l 1 -0 A + 


350905.1 


2.3 + 0.3 


-0.6 + 0.3 


33.1 + 0.8 


26.5 + 0.6 


75.68 + 0.07 


CH3OH 


7,, = 4! - 3 ^ 


358605.8 


1.0 + 0.3 


-2.1 +0.2 


24.5 + 0.8 


15.9 + 0.5 


0.36 + 0.05 


CH3OH 


7* = 7 2 - 6x £ 


363739.8 


1.7 + 0.3 


-0.4 + 0.3 


19.0 + 0.8 


15.1 +0.6 


143.75 + 0.13 


CN 


Nj = 3 7/2 - 2 5/2 AF=0 


340263.7* 


-1.3 + 0.4 


-5.4 + 0.3 


22.3 + 1.1 


17.3 + 0.8 


3.03 + 0.13 


CN 


Nj = 3 5/2 - 2 3/2 AF=1 


340033.5'' 


-1.4 + 0.1 


-4.8 + 0.1 


24.0 + 0.3 


19.4 + 0.2 


105.12 + 0.04 


CN 


tf, = 3 7/2 - 2 5/2 AF=1 


340248.2* 


0.2 + 0.2 


-5.7 + 0.1 


27.5 + 0.4 


20.7 + 0.3 


93.23 + 0.03 


cs 


7 = 7-6 


342882.9 


0.1 +0.1 


1.3 + 0.1 


31.3 + 0.4 


21.1 +0.2 


82.87 + 0.02 


c 34 s 


7 = 7-6 


337396.5 


-4.0 + 0.4 


0.7 + 0.3 


20.9 + 1.0 


13.0 + 0.6 


117.48 + 0.07 


o-H 2 CO 


Jk g ,k c = 5i,5 - 4i s 4 


351768.6 


-0.8 + 0.2 


-2.7 + 0.1 


28.1 +0.4 


19.9 + 0.3 


179.56 + 0.03 


p-H 2 CO 


Jk u ,k c - 5o,5 - 4 ,4 


362736.0 


-0.6 + 0.2 


-2.7 + 0.1 


22.8 + 0.4 


16.1 +0.3 


169.49 + 0.04 


p-H 2 CO 


Jk q ,k c = 5 2 ,3 - 4 2;2 


365363.4 


1.4 + 0.2 


-1.5 + 0.2 


21.1 +0.5 


14.0 + 0.3 


20.99 + 0.04 


HCN 


7 = 4-3 


354505.6' 


1.9 + 0.1 


-0.9 + 0.1 


26.5 + 0.3 


19.4 + 0.3 


95.18 + 0.03 


H 13 CN 


7 = 4-3 


345339.7 m 


2.4 + 0.1 


-3.3 + 0.1 


17.0 + 0.2 


14.0 + 0.2 


127.44 + 0.05 


HCO + 


7 = 4-3 


356734.2 


-1.2 + 0.1 


-1.5 + 0.1 


31.3 + 0.3 


22.2 + 0.2 


75.62 + 0.02 


H 13 CO + 


7 = 4-3 


346998.3 


0.6 + 0.1 


-0.8 + 0.1 


23.7 + 0.3 


18.2 + 0.3 


93.28 + 0.04 


HNC 


7 = 4-3 


362630.3 


-2.1 +0.1 


-4.1 +0.1 


25.8 + 0.2 


17.4 + 0.1 


75.90 + 0.01 


N 2 H + 


7 = 4-3 


372672.5 


-3.4 + 0.3 


-3.1 +0.2 


40.1 +0.9 


17.7 + 0.4 


69.61 + 0.02 


S0 2 


Jk u ,k c - 23 3 , 2 i - 23 2 , 22 


336089.2 


2.1 +0.4 


-2.1 +0.5 


18.8 + 1.2 


12.8 + 0.8 


55.95" +0.11 


so 2 


^ fl ,^ c ~ ^3,3 _ 4 2;2 


351257.2 


-0.3 + 0.1 


-3.1 +0.1 


18.0 + 0.2 


13.8 + 0.2 


141.41 +0.04 


so 2 


Jk u ,k c = 23 2 , 22 - 23i 5 23 


363925.8 


0.8 + 0.4 


0.9 + 0.5 


21.4+ 1.3 


10.5 + 0.6 


58.52" + 0.05 


so 


Nj = 7 B - 6 7 


340714.2 


3.1 +0.1 


-5.3 + 0.1 


17.8 + 0.3 


15.2 + 0.3 


118.64 + 0.08 


so 


= 8 8 - 7 7 


344310.6 


2.0 + 0.1 


-2.8 + 0.1 


18.7 + 0.2 


16.9 + 0.2 


152.92 + 0.07 


so 


Nj = 9 8 - 8 7 


346528.5 


1.3 + 0.1 


-1.7 + 0.1 


22.0 + 0.3 


17.5 + 0.2 


109.09 + 0.03 


continuum 


850 //m 


(N/A) 


-1.2 + 0.1 


-1.0 + 0.1 


36.6 + 0.2 


27.4 + 0.2 


65.51 +0.01 



Notes. Spatial extent determined from two-dimensional Gaussian fit to integrated intensity map of each transition. Uncertainties for peak position, 
major and minor axes, and position angle quoted here ar e formal fitting e rrors and do not include systematic errors due to pointing uncertainties. 
(a) Laboratory frequency of line transition from CDMS ( Miiller et al.|2005| ). 

{b) Offsets in right ascension (ARA) and declination (ADec) with respect to the reference position 20 h 29 m 24 s .9, 40 o ll , 19 ,, (J2000). 

(c) No beam deconvolution is applied before determining angular sizes. 

(d) p os iti on angle of the fitted major axis, expressed from north to east. 

(e) Blend of fourteen hyperfine components (total separation 1.6 MHz). 
(/) Blend of F=5-4 and 4-3 hyperfine components (separation 1.3 MHz). 
(g) Blend of F=4-3 and 3-2 hyperfine components (separation 1.4 MHz). 

Blended with 7_ 2 - 6_ 2 transition (separation 1.3 MHz). The frequency listed in the table is the average. 
(z) Blended with OS 18 0. Listed frequency is the average for the OS 18 5 3j 3-4 2;2 transition and the F=5/2-5/2 and 7/2-7/2 hyperfine components 
of CN 3 7/2 - 2 5/2 (total separation 3.2 MHz). 

0) Blend of F=7/2-5/2, 3/2-1/2 and 5/2-3/2 hyperfine components (total separation 4 MHz). 
{k) Blend of F=7/2-5/2, 9/2-7/2 and 5/2-3/2 hyperfine components (total separation 0.8 MHz). 
(Z) Blend of F=4-4, 3-2, 4-3, 5-4, 3-4 and 3-3 hyperfine components (total separation 3.6 MHz). 

(m) Blend of S0 2 13 2 ,i2 - 12 U i and hyperfine components F=4-4, 3-2, 4-3^ 5-4, 3-4 and F=3-3 of H 13 CN 4-3 (total separation 3.6 MHz). 
(n) Position angle poorly represent actual morphology of the map (see Fig.pl). 



N2H + shows the same on the northern side, giving rise to a four- 
directional morphology. 

The emission from HNC 4-3 and from the three CN tran- 
sitions (Fig. [2]) peaks at a position offset ~3" to the southwest 
with respect to the dust emission and the other molecular maps. 
Most other molecular maps have peak positions at ADec be- 
tween -3" and -1", and the 850/mi dust emission peaks at 
ADec^-rXTablefT]). Moreover, the peaks of the SO maps, espe- 
cially 7g-6 7 , appear to be shifted -4-5" to the southeast (Fig. [3}. 



In contrast, the CS and C 34 S maps (Fig. [3]) have peak positions 
offset 2-3" to the north. Finally, the two C2H maps, although 
not obvious from their fitted peak positions, are displaced to 
the northwest (Fig. [4]). We note that these position offsets are 
only marginally significant compared to the pointing accuracy 
(Sect. [2]). 

A spatially resolved two-peak structure is apparent in the 
map of CH3OH \2\-\2q (£ up =197 K) in the top right panel of 
Fig. [4] The two peaks are separated by roughly 12" (12 000 AU). 
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Table 2. Maximum integrated intensity values and velocity parameters. 



Molecule 


Transition 


J7 Ih 


integration 


maximum 


Gaussian line profile 








interval 61 


j T mh dV b 


^centroid 


FWHM J 






(K) 


(km s ) 


(Kkms -1 ) 


(kms -1 ) 


(kms -1 ) 


CO 


7 = 3-2 


33.2 


[-29, +6] 


410.9 


-4.36 + 0.29 


4.8 + 0.6 e 


13 co 


7 = 3-2 


31.7 


[-18, +6] 


225.5 


-5.92 + 0.11 


5.4 + 0.3 e 


c 17 o 


7 = 3-2 


32.4 


[-16, +4] 


18.6 


-5.86 + 0.06 


3.6 + 0.2 e 


C 2 H 


Nj = 4g/2 - 3 7 / 2 


41.9 


[-20, +6] 


17.2 


-5.99 + 0.09 


4.3 + 0.3 e 


C 2 H 


Nj = 4 7/2 - 3 5/2 


41.9 


[-15, +4] 


14.0 


-5.81 +0.13 


4.8 + 0.3 e 


CH3OH 


J K = 7i- 61 A + 


79.0 


[-12,-1] 


2.7 


-5.92 + 0.12 


4.0 + 0.3 


CH3OH 


7 i , = 12 1 -12 A T 


197.1 


[-11,0] 


1.7 


-5.36 + 0.19 


3.5 + 0.5 


CH3OH 


Jk = 7 - 6 E 


78.1 


[-10,-1] 


2.9 


-5.60 + 0.11 


3.5 + 0.3 


CH3OH 


J K = 1. 1 -6. 1 E 


70.5 


[-11,0] 


3.3 


-5.93 + 0.05 


3.3 + 0.2 


CH3OH 


Jk = 7 - 6 A + 


65.0 


[-14, +1] 


4.1 


-5.80 + 0.07 


3.7 + 0.2 


CH3OH 


J K = 1 2 -6 2 E 


87.3 


[-11, +1] 


3.1 


-5.60 + 0.07 


3.8 + 0.2 


CH3OH 


7„ = 4 - 3_i E 


36.3 


[-16, +1] 


5.8 


-7.01 +0.13 


5.7 + 0.3 


CH3OH 


J K = h-0 A + 


16.8 


[-11,0] 


3.5 


-6.05 + 0.05 


3.5 + 0.2 


CH3OH 


J K =4 Y -%E 


44.3 


[-11,0] 


3.5 


-6.42 + 0.07 


3.8 + 0.2 


CH3OH 


Jk = 1i-6iE 


87.3 


[-10,0] 


2.6 


-5.55 + 0.09 


3.2 + 0.3 


CN 


Nj = 3 7/2 - 2 5/2 AF=0 


32.7 


[-13,-1] 


3.3 


-5.38 + 0.24 


6.6 + 0.6 


CN 


Nj = 3 5/2 - 2 3/2 AF=1 


32.6 


[-16,0] 


21.5 


-5.80 + 0.19 


6.6 + 0.5 


CN 


Nj = 3 1/2 -2 5/2 AF=l 


32.7 


[-14, +2] 


25.7 


-5.61 + 0.08 


4.0 + 0.2 


cs 


7 = 7-6 


65.8 


[-14, +3] 


27.3 


-5.91 + 0.04 


3.5 + 0.1 e 


c 34 s 


7 = 7-6 


50.2 


[-12,-1] 


2.6 


-5.68 + 0.07 


3.3 + 0.2 


o-H 2 CO 


Jk u ,k c - 5i,5 - 4 1>4 


62.5 


[-13, +2] 


17.9 


-5.91 + 0.07 


3.9 + 0.2 


p-H 2 CO 


•^,,£c - 5o,5 - 4 ,4 


52.4 


[-11,0] 


8.3 


-5.62 + 0.05 


3.4 + 0.1 


p-H 2 CO 


•^r a ,^c ~ ^2,3 _ 4 2;2 


99.7 


[-13, +2] 


3.5 


-5.48 + 0.09 


3.8 + 0.3 


HCN 


7 = 4-3 


42.5 


[-17, +5] 


67.5 


-5.68 + 0.06 


5.3 + 0.2 e 


H 13 CN 


7 = 4-3 


41.4 


[-15, +3] 


10.5 


-5.32 + 0.07 


5.8 + 0.2 


HCO + 


7 = 4-3 


42.8 


[-18, +5] 


100.7 


-6.10 + 0.08 


5.2 + 0.2 e 


H 13 CO + 


7 = 4-3 


41.6 


[-13,0] 


7.9 


-6.03 + 0.07 


3.1 +0.2 


HNC 


7 = 4-3 


43.5 


[-12, +1] 


21.4 


-5.67 + 0.05 


3.4 + 0.1 


N 2 H + 


7 = 4-3 


44.7 


[-11,-2] 


9.7 


-5.66 + 0.06 


3.1 +0.2 e 


S0 2 


Jk u ,k c - 233 ; 2i - 23 2j22 


276.0 


[-11,-1] 


2.6 


-5.64 + 0.17 


3.2 + 0.4 


so 2 


Jk ,k c - ^3,3 - 4 2;2 


35.9 


[-11,0] 


4.5 


-5.65 + 0.06 


5.3 + 0.2 


so 2 


Jk u ,k c - 23 2j22 - 23i ? 23 


259.9 


[-11, +1] 


2.2 


-5.37 + 0.17 


4.7 + 0.4 


so 


Nj = 7 8 - 6 7 


81.2 


[-13, +2] 


11.9 


-6.01 + 0.05 


5.2 + 0.2 


so 


= 8 8 - 7 7 


87.5 


[-14, +1] 


11.0 


-6.08 + 0.05 


5.2 + 0.2 


so 


Nj = 9 8 - 8 7 


78.8 


[-15, +5] 


14.0 


-5.91 +0.17 


5.8 + 0.4 



Integration is done including the entire channel in which the boundary falls. 

Uncertainties in absolute intensity calibration are 15%. Maximum intensity levels are used to scale the contour levels in Figs.[2}j4] 
Centroid velocity, with respect to Vi sr , determined from Gaussian fit to the velocity profile at the central spatial position. 
Full width at half maximum determined from a Gaussian fit to the velocity profile at the central spatial position. 
Apparent non-Gaussian line profile and/or presence of line wings outside of Gaussian. 



The southeast warm CH3OH peak is associated to the infrared 
source detected by | Van der Tak et aL ( |1999| their Fig. 1), the 
canonical center of the envelope o f dus t and gas. We discuss this 
two-peak structure further in Sect. |5.2| The only other molecular 
transitions shown in this paper with an upper level energy above 
100 K are the two SO2 23-23 transitions, where a double peak 
structure is not seen. The N2H + morphology does not coincide 
with the spatially separated warm methanol pockets, neither in 
angular separation, nor in position angle. 



The high-excitation 7=23-23 transitions SO2 transitions in 
Fig. [3] appear to be stretched in a direction perpendicular to the 
outflow direction. The fitted position angle for these two maps 
does not reflect the actual orientation of the extended emission 
(cf. Table [T] and Fig. [3]), perhaps due to the boxy non-Gaussian 
shape of the emission. The real position angle for these maps 
would plausibly be around 140°-150°. 



3.3. Velocity structure 

The velocity structure of the molecular gas can be illustrated 
by line profiles in velocity space, and by 'channel maps'. First, 
the velocity profiles of all 35 lines at the central spatial position 
are fitted by a Gaussian to determine the centroid velocity and 
linewidth (Table [2]). We find centroid velocities between -5 and 
-7kms _1 , with a median of -5.8 km s" 1 . This is slight ly lower 
than the previo usly determined velocity of -5.5 km s" 1 (Van der 
Tak et al.|l999| . Gaussian linewidths range from 3 to 7 kms _1 , 
with a median of 3.9 km s" 1 . 

Second, to show the velocity structure of the overall molec- 
ular envelope, each integration interval from Table [2] is di- 
vided into three velocity bins: a central bin defined as the 
[-1.5, +1.5] km s" 1 interval, relative to the systemic Vi sr of 
-5.8 km s" 1 , and 'blue' and 'red' bins composed of all emis- 
sion below and above the central velocity bin, respectively. The 
resulting channel maps are presented in Figs. [6] and [7] We only 
show a map of velocity structure for those transitions that show 
significant emission in at least two of the three velocity bins. 
In addition, we present a velocity map where emission is di- 
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28s 24s 20h29m20s28s 24s 20h29m20s28s 24s 20h29m20s28s 24s 20h29m20s 

RAQ2000) RAQ2000) RA (J2000) RA Q2000) 

Fig. 6. Channel maps for the CO isotopologs, N-bearing species, HCO + and H 2 CO. The velocity range is divided in three ve- 
locity bins with respect to the systemic velocity of -5.8 km s" 1 : [-1.5, +1.5] km s" 1 is shown in grayscale stretched between 
0.7 K km s" 1 (white) and the maximum intensity in the central channel (black); all emission at relative velocities below -1.5 km s" 1 
is shown in blue contours, and all emission at relative velocities above +1.5 km s" 1 in red contours. The lowest contour level is 
manually fine-tuned to emphasize velocity structure in each case without showing too much noise. Consecutive contour levels are 
drawn in steps of 10% of the maximum intensity in the brightest channel. Lowest contours: 15.0 K km s" 1 for CO, 5.0 K km s" 1 for 
13 CO, O^Kkms" 1 forC 17 0, 0.35 K km s" 1 for HCN and HNC, O^Kkms" 1 for N 2 H + , 0.8 Kkms" 1 forCN 3 5/2 -2 3/ 2, l.OKkms" 1 
for CN 37/2-25/2, and 0.8, 0.4 and 0.3Kkms _1 for the three respective H2CO transitions. The FWHM of the telescope beam is 
indicated in the top left corner of each panel. 



vided over seven velocity bins for the bright 13 CO 3-2 transition 
(Fig.[8]). 

A general west-east gradient from blueshifted (approaching) 
to redshifted (receding) emission is perceived in the CO species, 
N 2 H + , CN, HCO + , HCN, and HNC (Fig. [6]). The same trend 
is seen in the red and blue peak positions in maps of C2H and 
CS (Fig. [7]). Although the typical angular separation between 
the peak of the blue and red components at <10" is less than 
the FWHM of the telescope beam, the high signal-to-noise ratio 
in our maps allows for a relative position determination down 
to a fraction of the beam size. Moreover, the total extent of the 
blue and red lobes for the molecules mentioned above is signif- 
icantly separated, with red contours reaching positions > 30" 
east, and blue contours > 30" west of the central position. 
Conversely, blue and red components are coincident for SO and 
SO2. H2CO maps show no clear velocity pattern. Velocity struc- 
ture in maps of CH 3 OH is complicated by blueshifted emission 
from the northeastern plume (see below, and Sect. |3.2| ). Velocity 
structure in the spatially separated peaks of CH3OH 12i-12q is 



not discernable in our observations. Hence, the line of sight ve- 
locity of the two pockets must be within 1 km s" 1 of each other. 
In addition, the majority of the velocity maps exhibits emission 
that is stronger in the blueshifted bin than in the redshifted bin. 
Exceptions are HNC 4-3 and the P-H2CO transitions, where the 
red emission is stronger than the blue emission, and the SO tran- 
sitions, where intensity levels in the red and blue bin are compa- 
rable. 

Upon inspection of the Gaussian line profiles parametrized 
by Vcentroid and FWHM in Table [2j the full integration interval is 
notably wider than the Gaussian profile for the CO species, as 
well as for C 2 H, CS, HCN and HCO + . This is due to our orig- 
inal choice to incorporate all emission from a particular molec- 



ular transition in the corresponding J T m ^dV map (Sect. 3.1). 
The non-Gaussian part of the line intensity at the central spaTTal 
position mainly lies on the blueshifted side for these lines (see 
also Fig. [9] and the notes in Table [2]). For most of the lines men- 
tioned above the line wings are relatively weak, comprising up to 
20% of the total integrated line intensity, but the single Gaussian 
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+40 




+40 



28s 24s 20h29m20s 

RAQ2000) 



+40°10'20" 



28s 24s 20h29m20s28s 24s 20h29m20s28s 24s 20h29m20s 

RAQ2000) RAQ2000) RA Q2000) 



Fig. 7. Channel maps for sulfur bearing species, CH 3 OH, and C2H. Colors and contours have the same meaning as in Fig.pJ except 
lowest contours: O^Kkms" 1 for SO 7 8 -6 7 , O^Kkms" 1 for SO 8 8 -8 7 and 9 8 -8 7 , O^Kkms" 1 for CS 7-6, 0.2Kkms" H for S0 2 
53 9 3-4 2 ,2, 0.35, 0.6, 0.35 and 0.3 K km s" 1 for the four consecutive CH 3 OH transitions, and 0.4 K km s" 1 for both C 2 H transitions. 



profile is an especially poor representation for the main isotopic 
CO 3-2 line profile, where as much as 45% of the integrated line 
intensity falls outside of the fitted Gaussian. We note that wings 
might be present even in other (weaker) lines, which would be 
invisible in the noise of the observations (Sect. [2]). 

The velocity components of the N 2 H + 4-3 line are special in 
the sense that they are spatially resolved: the emission peaks of 
the blue and red N 2 H + wings in Fig.[6]are separated by more than 
one beam size. For other species, the different components over- 
lap in velocity space as well as spatially. Like for several other 
species, emission is significantly stronger in the blueshifted ap- 
proaching component than in the recedi ng c omponent. In the 
four- arm structure noted for N 2 H + in Sect. |3.2| the southeast arm 
consists exclusively of emission from the red velocity bin. 

The northern plume noted in Sect. 3.2 has more contribution 
from the blueshifted component of the emission than from the 
redshifted and central velocity bins, as shown in Figs. [6] and [7] 
This is particularly clear in velocity maps of CH3OH. 

4. Radiative transfer modeling 



Section |4.1| deals with dust continuum modeling and describes 
the derived density structure of the envelope. The spatial distri- 
bution of molecular emission from a radiative transfer model is 
addressed in Sect. A2_ Section [43] explores the effects of differ- 



ent density structures on molecular emission. Finally, the model 
from Sect.P 



4.2 



ture (Sect. 



is adjusted to examine the effects of velocity struc- 



4.4) and a flattened morphology (Sect. 03 



4.1. Dust continuum models 



We use the dusty code ( flvezic et al.||1999| ) to create a grid of 
dust continuum models starting from the density profile of mate- 
rial in a spherical (n oc r~ a ) envelope around a luminous source. 
Parameters that vary across the grid include the temperature at 
the inner edge of the dust envelope, and the power law index 
a of the density distribution. The internal luminosity is fixed at 
2xlO 4 L at lkpc distance. A self-consistent dust temperature 
and radiation field is derived within the envelope, taking into 
account scattering, absorption and emission processes along a 
radial path. The intrinsically one-dimensional models are resam- 
pled onto a radial grid of sky coordinates to obtain an intensity 
map, which is subsequently convolved with a synthetic telescope 
beam. 

The resulting radial intensity profile of each model is com- 
pared with the 450 and 850 /mi SCUBA continuum maps 
(Sect. [2]). The best-fit dusty model for 450 yum has a = 1.10, 
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28s 24s 20h29m20s28s 24s 20h29m20s28s 24s 20h29m20s 

RAQ2000) RAQ2000) RA Q2000) 

Fig. 8. Emission from the 13 CO 3-2 transition divided over seven velocity bins. Line intensity is summed over the velocity range 
indicated in each panel (in km s" 1 ) and represented as black contour lines. Contour levels are 90%, 80%, . . . , 10% of the maximum 
integrated intensity in the velocity bin: 3.9, 14, 41, 124, 50, 8.1, and 0.95 Kkms" 1 in the respective bins. Like in Fig. [2] the 50%- 
level contour is thicker, and any contours below 0.7 Kkms -1 are not drawn. The total integrated 13 CO 3-2 line intensity over the 
entire velocity range is shown in colorscale in each panel, for reference. 



whereas the one for 850 yum has a = 1.45. We take the modeled 
dust opacity, r y , and use the relation 



,p(r)&r 



(1) 



to derive the scaling factor p in the mass density, p(r) = 
Po( r I ' r o)~ a • Here, we assume values for the absorption coeffi- 
cients #850 = 0.03 cm 2 g" 1 and K450 = 0.07 cm 2 g" 1 . The mass 
density scaling factor p is divided by the mass of a hydrogen 
molecule in order to convert to gas number density (no). The de- 
scription ofdensityand temperature used for the a = 1.5 models 
in Sects. [431 and |4.4| is derived from an a - 1.50 dusty model, 
which is very similar to the best- fit solution with a - 1.45. 



I Van der Tak et al.| ( |2000| ) constrained the slope of the power 
law density profile, a, using CS line observations. With the index 
being allowed to take on values of 0.5, 1.0, 1.5 or 2.0, they found 
that the best fit to several observed CS lines is given by a = 1.0. 
The observed 450 /mi continuum map also shows good corre- 
spondence to the radial intensity profile of an a = 1.0 model, but 
the 850 fim co ntinuum map is fit better with a = 1.5 (Sect. 5.3 of 
Van der Tak et al.|2000 ). We have independently confirmed the 
result that 850 fim dust emission favors a steeper density gradi- 
ent^ = 1.45) than 450 yum (a = 1.10). In addition, [De Wit et al.| 
( 2009 ) find that a power law density profile with a = 1.0 fits 
the radial intensity profile of high-resolution (0.3" half-power 
beamwidth) 24.5 //m observations, but they also report that a fit 
to the spectral energy distribution prefers a = 1.5. Based on the 
above arguments, we choose a density profile with a = 1.0 for 
the molecular radiative transfer in Sects. |4~2]|4.4 4.5 An explo- 
ration of different values of a is addressed in Sect. 14.31 



4.2. Static spherical model 

In the quantitative modeling of the physical structure of the en- 
velope that gives rise to the molecular emission, we focus on 
C 17 0, HCN, CS, C 34 S, HCO + and H 13 CO + . This set of species 
is chosen to cover a variety of atomic constituents, and to have 
observed line strengths that provide the highest signal-to-noise 
available in our sample. While CS and HCO + are juxtaposed 
with rarer isotopic variants, the isotopolog H 13 CN is not in- 
cluded because its line signal is blended with an SO2 transi- 
tion (see Table 0). The maps of C 17 0, HCN, CS and HCO + 
have observed major to minor axis ratios below 1.5 (Table [TJ, 
and they have a largely homogeneous spatial distribution. This 
makes them suitable t racer s of the quiescent overall envelope 
of AFGL2591. Section 5.2 addresses the interpretation of distri- 
butions for species that are not explicitly modeled here but do 
appear in Table [T] 



We use the radiative transfer code ratran (Hogerhei jde & 
|Van der Tak 2000) to test one-dimensional source models for 
AFGL2591. This relies on molecular data files from lamda 



([Schoier et al.|2005 ), with collisional rate coefficients calculated 



by|Green & Th addeus ( 19741 HCN),|Turner et al.| ( [T992| CS), 
Flower| ( |1999| HCO + ), and |Yang et al.|(|2010| CO). Motivated 



by the dust modeling described in Sect. |4.1| we adopt a physical 
model from I Van der Tak et al.| fL999{ |2000), which assumes a 
power law density profile of the form: 



(2) 



n(r) = n — 
^0 



with ft =5.3xl0 4 cm -3 at r =2.7xl0 4 AU, and a=l. The model 
cloud is truncated at a radius of 6.2xl0 4 AU (62"), and has a to- 
tal H2 mass of 193 M©. The innermost shell of our model extends 
from r = to -200 AU, and has an H 2 density of 1.2xl0 7 cm" 3 



10 



M. H. D. van der Wiel et al.: JCMT-SLS: The molecular envelope of AFGL2591 



and a temperature of 372 K. The density singularity at r = 
is thus avoided by the discrete sampling over our finite number 
of model shells. A centrally peaked temperature profile ranging 



down to 25 K was derived from dust emission by | Van der Tak 



et al. (2000). Molecular abundances are kept constant as a func- 



tion of position in the envelope. A constant turbulence is incor- 
porated, parametrized by the 1/e halfwidth velocity of 2kms _1 , 
consistent with the typical value of 3-4 km s" 1 of the measured 
FWHM linewidths (Table[2]). Initially, no infall, outflow, rotation 
or other velocity gradients are included. We refer to this model 
as the 'static spherical model'. 

We use the ray tracing routine of ratran to produce spectral 
maps with 0.5" pixels at the distance of 1 kpc. These maps are 
intentionally oversampled with respect to the spatial resolution 
of our observations, and they are subsequently convolved with 
a Gaussian beam profile appropriate for the molecular transition 
under consideration. 

The spatial distribution of observed integrated line emission 
is probed by a slice along the 'long axis' (position angle 75°) and 
another slice along the 'short axis' (position angle 165°) of the 
flattened large scale cloud. These slices are compared with the 
static spherical model in Fig. 10 The two observed slices pro- 
vide a measure of the non- spherical morphology of the cloud in 
the plane of the sky. Any modeled intensity profile which differs 
from the observed profiles by more than the difference between 
the observed 'long axis' and 'short axis' profiles should be re- 
garded as incorrect. 

We attempt to match the modeled peak line strengths to 
the observed ones by adjusting individual molecular abun- 
dance values, keeping isotopic ratios fixed at CS/C 34 S=20 and 
HCO + /H 13 CO + =60. We find abundances, relative to H 2 , of 
2xl0" 9 for C 34 S and 2.5xl0" 10 for H 13 CO + . Consequently, the 
abundances of CS and HCO + are then fixed at 4xl0" 8 and 
1.5xl0~ 8 , respectively. Since we do not have an unblended op- 
tically thin isotopolog for HCN, we take 2xl0" 8 ; for C 17 we 
adopt 5xl0~ 8 as an abundance. These abundance values a re con- 
sistent with earlier estimates by Van der Tak et al. ( 1999 ). 

The shapes of the modeled intensity profiles as a function of 
radius correspond to the observations as long as the line emission 

^C 34 S,H 13 CO + ),butthe 
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is predominantly optically thin (Fig. ^ 
emission of both HCN and HCO + grows optically thick before 
the modeled line intensity reaches the observed value. This re- 
sults in modeled emission profiles which are too weak near the 
center of the cloud, and fall off less steeply than the observed in- 
tensity profile. A deviation from the fixed isotopic ratio is briefly 
investigated. However, since the cloud is already optically thick 
at the center, a higher molecular abundance for the main isotopic 
species does not increase the emerging line intensity at the cen- 
tral position, while the intensity at more optically thin lines of 
sight away from the center does increase, thereby adding to the 
mismatch of the overall spatial profile. Moreover, with observed 
variations of the isotopic ratio restricted to a factor < 3 (e.g., 
|Chin et al.|1996| ), our results are not affected. 

The modeled intensity of CS 7-6 is obtained by first adjust- 
ing the C 34 S abundance to match the observed peak intensity, 
and simply scaling the CS abundance by the canonical factor 20 
(see above). The modeled main isotopic CS 7-6 emission ap- 
pears to be stronger than observed, but we note that this discrep- 
ancy is well within the intensity calibration uncertainty of 15% 
of the observations. 

In the spectral dimension, each observed line has a roughly 
symmetric and single-peaked shape, as illustrated in Fig. [5] for 
HCO + 4-3. Spectral line shapes are similar for other species. 
Considering the line profiles resulting from the radiative transfer 



ID static model 

— infall model 
flat model at i=20° 




-20 -15 -10 -5 
y lsr (kms- 1 ) 



Fig. 9. Top: Velocity profiles of observed HCO + 4-3 (his- 
togram), compared to those from various models, all at the 
central spatial position. The profiles from the spherical (ID) 
model and the infall model follow from an HCO + abundance 
of 1.5xl0~ 8 , the abundance for the model used for the flattened 
model profile is 6xl0~ 9 . Bottom: velocity profiles of observed 
and modeled H 13 CO + , with abundances scaled by a factor 60 
with respect toH 12 CO + . 



models, individual optically thin lines (C 34 S 7-6, H 13 CO + 4-3, 
C 17 3-2) have very comparable shapes, matching observed line 
shapes. However, for the main isotopologs HCN and HCO + , the 
high optical depth at line center is prominently displayed in the 
form of a self-absorbed line profile (Fig. [9]). 

A particular explanation is required for the spatial distribu- 
tion of C 17 0, which is also discrepant w.r.t. the observed distri- 
bution but is optically thin in the modeled line. See Sect. |5.1| for 
a discussion. 

To help the line emission escape the molecular cloud, 
the opacity at line center might be diminished for HCN and 
HCO + by i ntrod ucing a velocity gradient in the protostellar 
cloud (Sect. 4.4 ), or by adjusting the geometry of the envelope 
(Sect. |4.5| ). In the originally optically thick cases this might re- 
sult in extra emerging line intensity. On the other hand, the in- 
tegrated line intensity of the optically thin H 13 CO + and C 34 S, 
for which the static spherical models already match the observa- 
tions, should be conserved. 
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Fig. 10. Position dependence of inte grate d line intensities for 
static spherical ratran models (Sect. |4.2) of C 17 0, HCN, CS, 
C 34 S, HCO + and H 13 CO + . Black solid and dashed lines repre- 
sent two perpendicular slice directions in the observed maps, red 
lines represent radial dependence of integrated line strength in 
the ID static models. At the line center, the sky ray tracing rou- 
tine calculates the following optical depths through the center of 
the model sphere: 0.095 for C 17 0, 33 for HCN, 10 for CS, 0.38 
for C 34 S, 34 for HCO + and 0.64 for H 13 CO + . Abundance values 
(X) are presented in the form a (b), representing ax 10^. 



4.3. Steeper density profiles 



We recall that we adopt a = 1 in Sect. |4.2| for the density pro- 
file given by Eq. ^ for the static spherical model. The dynamic 
spherical models (Sect. |4.4| ) and the flattened models (Sect. [43]) 
are based on the same density p ower law. The validity of this 
assumption is addressed in Sect. |4.1[ below we illustrate how 
steeper density profiles affect our results. 

Because the value of a is not unambiguously equal to 1.0, we 
run test models with a = 2.0 and a = 1.5 for the static spherical 
model for HCO + . Both models have self-consistent temperature 
profiles derived in earlier work ( [Van der Tak et al.||2Q00j > or by 



the dust modeling described in Sect. 4.1 The qualitative trend 



is that a steeper density profile results in a steeper temperature 
gradient. 

In the static a = 2.0 model, the HCO + 4-3 integrated line 
intensity emanating from the center is about 20% higher than 
that of the a = 1.0 model at the same abundance (cf. Fig. 10). 
However, it still only reproduces 60% of the observed peak inte- 
grated intensity. Doubling the HCO + abundance from 1.5xl0~ 8 



to 3x1 0~ 8 increases the central J T m ^dV value by only a few 
percent, with the optical depth at the model center at line cen- 
ter increasing from 28 to 52. Therefore, even with further in- 
creased abundance, the a = 2.0 model will still underproduce 
the observed peak J r m bdV level. On the other hand, the spatial 

profile of J r m bdV does fall off steeper than for a = 1.0, bring- 
ing the shape closer to the observations, but we note that the 
H 13 CO + J T mh dV profile becomes too sharp to match the obser- 
vations. Moreover, the fixed 1/60 ratio between the abundances 
of H 13 CO + and HCO + results in an overproduction of the cen- 
tral H 13 CO + J T m \>dV, simultaneously with the wftderproduction 
of HCO + emission. Thus, we conclude that a density distribu- 
tion with a = 2.0 does not improve the match to the observed 
intensity distribution with respect to the a = 1.0 case (Sect. 4.2 ). 

In the static a = 1.5 model, the HCO + 4-3 integrated line 
intensity d istrib ution is more peaked than in the a = 1.0 case 
from Sect. |4.2| However, the peak j r m bdV value, which was a 
factor ~2 too"low for a = 1.0, is increased by only <10% for 
a - 1.5. At the same time, H 13 CO + emission toward the central 
position is overproduced by a factor of 2. We conclude, as for 
a - 2.0 and 1.0, that a density power law slope of 1.5 fails to 
explain the spatial distribution of line emission if we consider an 
optically thick and an optically thin tracer simultaneously. 

4.4. Dynamic spherical models 

In this section we attempt to diminish the optical depth near the 
emission line frequencies by introducing velocity structure in the 
modeled cloud. Velocity structure can be considered as either 
random (turbulent) motion or systematic motion. We investigate 
the effects of a radial turbulence gradient and of a systematic 
infall signature on the line opacity and on the distribution of 
emerging intensity. 

While the originally constant turbulence of 2kms _1 is con- 
sistent with non-thermal line widths found in other massive star- 



forming regions ( Caselli & M yers 1995 ), we introduce a gradient 
ranging from lO krns" 1 in the central shell of the model down 
to lkms -1 in the outer shell. Even with this deliberately ex- 
treme turbulence gradient, the resultant line optical depth is not 
reduced significantly, whereas the line shapes become broader 
in velocity space, as expected. Since this broadening is inconsis- 
tent with the observed FWHM of ~4kms _1 for most lines, we 
refrain from further exploring this particular type of differential 
dynamics. 

Apart from a turbulence gradient, a systematic radial ve- 
locity gradient could be viable. If our envelope is to be form- 
ing a centrally condensed object, an inf ailing rather than an ex- 
panding velocity field seems reasonable. An infall signature is 
classically regarded as the double-peaked spectral line profile 
with a strong er blue than red peak in optically thick lines (see, 
e.g., Fig. 5 of |Evans|1999| ). Such a signature is either not present 
in our source, or is unresolved by the observations (spectral res- 
olution 1 kms -1 ). Any unresolved double-peak structure would 
have to be separated by <2kms _1 , yielding an upper limit to 
the infall veloci ty of 1 km s" 1 near the dust sublimation radius, 
-100 AU (e.g., |De Wit et al.||2009] ) ? close to the protostar. This 
limit is conservative, since even observations at slightly higher 
spectral resolution ( [Mitchell et al.|1992[|Van der Tak et al.|1999| ) 
show no sign of infall in line profiles of, e.g., 1J CO and CS. 
Choi et al.| ( [T995] ) and |Hogerheijde & SandeTTl ( [2000] ) have 



considered one or more off-center positions to constrain models 
of inf ailing protostellar envelopes. They find that inside-out in- 
fall models with a = 1.5 are among the possible fits for some 
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protostellar envelopes. We emphasize that different a values are 
found for different protostellar envelopes, although the power 
law index of the density profile a is generally between 1.0 and 
2.0. |Hogerheijde & Sandell| ([2000) highlight that the density and 
velocity structure are degenerate. 

To enable a direct comparison between results from the static 
and the inf ailin g mo dels, we retain the spherical a = 1.0 enve- 
lope from Sect. A2_ and add a radial velocity profile applicable 
for the freefall region of an inside-out collapsing cloud (follow- 
ing [Shu[f977]): 



VWallO) = - 



2GM(r) 



(3) 



where r represents radial distance from the cloud center and G 
is the gravitational constant. The enclosed mass M in Eq. ^ 
is defined as M(r) = M centra i + M env (< r). We set M centra i = 
16 M Q to account for the central protostellar object; M env (< r) 
represents the gas mass in the envelope contained in a sphere of 
radius r. The turbulent velocity is again constant at 2kms _1 in 
the infalling mo del. Relative molecular abundances are kept the 
same as in Sect. 14.21 

The velocity gradient has the desired effect of decreasing 
the optical depth at the line center and through the center of 
the spheric al m odel by a factor > 3. As in the spherical static 
case (Sect. |4.2| ), the modeled profiles of the rarer isotopologs 
C 34 S and H iJ CO + still match the observations after the intro- 
duction of infall. The spatial profile of J T m bdV of HCN 4-3 and 
HCO + 4-3 benefit from the reduced optical depth. However, it 
is insufficient to reproduce the observed peak J r m bdV value for 
HCN, as well as the spatial profile of HCO + . On the other hand, 
the infall model for CS 7-6 overproduces line intensity at the 
central position by -70%. In general, a significant discrepancy 
remains for the optically thick lines, both in terms of shape of the 
spatial profile as well as the peak integrated line intensity value. 
More importantly, the velocity gradient introduced in Eq. ([3} re- 
sults in a very asymmetric and double-peaked velocity profile 
for optically thick lines, and a wider profile (FWHM > 6 km s" 1 , 
Fig. [9]) for optically thin lines such as H 13 CO + 4-3. Concluding, 
the observed spectral line profiles indicate that an infall model is 
inadequate for the molecular envelope studied here. 

Finally, we test an a = 1.5 model for the infalling case, 
dynamically consistent with the physical description by Shu 



11) for this 



{ [1977) . The spatial distribution of f T mb dV (Fig. 
model comes closer than any other to a simultaneous - explana- 
tion of an optically thick (HCO + ) and an optically thin (H 13 CO + ) 
tracer, with abundances halved with respect to the a = 1.0 mod- 
els. However, when the spectral direction is also taken into ac- 
count, the same discrepancies emerge as for the other infall mod- 
els. We conclude that the adjustment of the density power law 
slope from 1.0 to 1.5 does not alleviate the mismatch between 
infall models and the observed spectral line profiles. In fact, the 
static models in Sects. |4.2| and |4.5| produce better matches to the 
observed spectral line profiles than the infall models. 

4.5. Static model with flattened geometry 

After considering a static sph erica l model (Sect. |4.2| ) and dy- 
namic spherical models (Sect. 4.4 ), we now explore a morpho- 
logical model of the molecular envelope which is flattened in 
one direction. The spatial axes of the model perpendicular to 
the line of sight are constrained by the observed morphology of 
the molecular cloud, which indicates a projected cloud structure 
not far from circular (major to minor axis ratios < 1.5). On the 
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Fig. 11. Position dependence of integrated line intensities for 
HCO + an d H 1 3 CO + 4-3 resulting from an infall model with a = 
1.5 (Sect. [43] ). Optical depths at line center through the center 
of the model cloud are 4.6 for HCO + and 0.06 for H 13 CO + . 



other hand, the cloud can be flattened significantly in the direc- 
tion along the line of sight. Moreover, if a small inclination of the 
short axis of the modeled system w.r.t. the line of sight ('view- 
ing angle') is introduced, we expect to be able to explain part of 
the non-circular geometry in the observed maps. Geometrically, 
the observed major/minor axis ratios on the sky of < 1.5 would 
constrain a flattened structure to be positioned at an inclination 
angle < 45°. 

For the density structure of the flattened model, we take an 
ellipsoidal adaptation of the spherical description from Eq. ([2]), 
where two of the three spatial axes are identical to the spheri- 
cal case and the third axis is compressed. The spherical density 
structure is first translated into cylindrical coordinates (R, z) and 
subsequently flattened by compressing the z-direction by a factor 
of 3. The density scaling factor no is scaled up by the same fac- 
tor, such that the total molecular mass of the model is conserved. 
The temperature structure in the cloud is taken to depend on the 
cumulative column density toward the central p ositio n, which 
is extracted from the one-dimensional case (Sect. |4.2| ). There is 
no bulk velocity structure in this model. The turbulent veloc- 
ity is constant at 2kms _1 , as in the s phe rical static and spher- 
ical infalling models from Sects. 4.2 and 4.4 For the radiative 
transfer, we use the axisymmetric version of the ratran code 
( [Hogerheijde & Van der Tak|2000| ). 

Figure [12] presents the result of radiative transfer using the 
ellipsoidally flattened model for HCO + 4-3. We choose two rel- 
ative molecular abundances, X = 1.2xl0~ 8 and 6.0x1 0~ 9 , and 
three inclination angles for the ray tracing, i = 20,40, and 60°, 
where 0° is 'face-on' (viewing along the z-axis) and 90° is 'edge 



on'. The best match is given by X = 6xl0 -9 , i = 20° (Fig. [12 
bottom right panel). The figure shows that in all cases, the mod- 
eled HCO + 4-3 distribution is much more extended than what 
is observed. We note that the flattened envelope at X = 6xl0~ 9 
and i = 20° does explain the observed J T m \>dV of HCO + 4-3 
at the central position, with a molecular abundance of 6xl0~ 9 , 
whereas the spherical models (Sect. A2_ AA\ underpredict the 
observed J T m \>dV at the central position, even with more than 
double the abundance (X = 1.5xl0~ 8 ). The flattened model also 
has a significantly lower optical depth of 13 (X=6xl0~ 9 , /=20°) 
at line center through the center of the model, where the spheri- 
cal models have r > 20 for HCO + 4-3. For i = 40°, at the same 
low abundance, the optical depth is 15, and for i = 60° it is 23, 
comparable to opacities in the spherical models. 

The modeled HCO + 4-3 line shape in velocity space is also 
discrepant w.r.t. the observed line profile, which has a roughly 
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Fig. 12. Position dependence of integrated line intensities for 
HCO + 4-3, compared with integrated line intensity profiles from 
flattened ellipsoidal physical models (Sect. |4.5b . Model slices 
are shown for an HCO + abundance of 1.2x10"^ (left column) 
and 6xl0 -9 (right column). The viewing angle i varies from 60° 
(top) to 40° (middle) to 20° (bottom). Values of optical depth 
at line center through the center of the model, as determined 
by the ray tracing routine, are for X - 1.2xl0 -8 : 42, 28, 23 at 
i = 60,40,20°, respectively; and f or X = 6xl0" 9 : 23, 15, 13 at 
i = 60, 40, 20°, respectively. 



Gaussian shape with one single peak centered at the known 
source velocity. Although an optical depth of 13 at line center 
is lower than in the spherical models, Fig. [9] shows that the flat- 
tened model still results in a very thick line profile, where most 
of the photons find a way out of the cloud at > 2kms _1 away 
from the line center. As seen before, the modeled 4-3 line of the 
isotopic variant H 13 CO + is optically thin (r = 0.3), and there- 
fore produces a better match to the observed line profile (Fig.[9j 
bottom panel). 

To test how sensitive the emerging integrated intensity pro- 
files are to the adopted flattening factor of 3, we also investigate 
the J 7" m bdV map of HCO + 4-3 resulting from a model with a 
flattening factor of 6. The peak observed J T mh dV value is now 
best matched with an HCO + abundance of ~2xl0" 9 (at i = 20°), 
with an optical depth at line center of ~5. As previously, the 
spatial profile of J r m bdV is too shallow with respect to the ob- 
served map. Hence, model geometries with flattening factors of 
3 and 6 both fail to explain the observed intensity distribution 
of HCO + 4-3. While we consider flattening factors of a few to 
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Fig. 13. Relative level populations of the / = 3 and J = 2 levels 
of C 17 0, and the / = 4 level of H 13 CO + , as determined by the 
Monte Carlo method in ratran for the static spherical model. 
The population in each / level is given relative to the total pop- 
ulation for that species. The gas temperature profile (T^ n ) of the 
model cloud is shown by the solid blue line. The energy level 
of the / = 3 rotational state of C 17 is indicated by the dotted 
horizontal line. 

be reasonable, much flatter models would result in a geometrical 
description of the envelope that resembles a sheet, for which no 
supporting evidence exists. 

5. Discussion 

5.1. Physical structure of the modeled envelope 

In Sect. [4] we have compared the observed spatial distribution of 
a selection of molecular tracers with spherical static and infalling 
models, and with a static ellipsoidally flattened model. The spec- 
tral lines of the rare isotopic variants C 34 S and H 13 CO + can be 



explained by the static spherical model from Sect. |4.2| both in 
terms of integrated intensity at emission center and shape of the 
spatial profile. But for the main isotopic species HCN, HCO + , 
and CS, as well as for C 17 0, this representation of the molecu- 
lar envelope fails to reproduce the observed integrated intensity 
levels. We conclude that for HCN, HCO + and CS the line optical 
depths of > 10 (see Fig. 10) are preventing radiation from escap- 
ing the model envelope. Attempts to compensate for the missing 
intensity by increasing the molecular abundance logically result 
in even higher line opacities and do not result in more photons 
escaping the cloud. Moreover, the already self-absorbed spectral 
line profiles become more discrepant w.r.t. the observations as 
the abundance is increased. 

Unlike the other two lines that show a large discrepancy in 
spatial distribution of the emission (HCN 4-3 and HCO + 4-3), 
the modeled C 17 3-2 emission is optically thin (r = 0.1) at 
line center through the center of the model cloud. The low opti- 
cal depth is confirmed by inspection of modeled line profiles in 
velocity space, which are well represented by a Gaussian, both 
on and off center. We suspect that the relative overproduction of 
C 17 3-2 away from the central position is caused by excitation 
effects. Figure 13 shows that the level populations of the / = 3 
level determined by the Monte Carlo method peaks at 1 8 000 AU 
(outside the FWHM of the central beam of 15" = 15 000 AU), 
and declines steeply toward smaller radii. If higher /-levels are 
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examined, it becomes apparent that the high kinetic temperatures 
in the center of the envelope leaves lower levels such as / = 3, 
with an energy of only 32 K above ground, sparsely populated. 

The mismatch with the observed spatial distribution of C 17 
could imply that either the description of the excitation condi- 
tions (i.e., temperature) is inadequate, or an additional compo- 
nent of cold, tenuous gas should be added to the model. The 
current truncation radius of 64kAU, at which point the temper- 
ature is 22.6 K, is derived from CS line observations by | Van der| 
Tak et al. ( 2000| ). An extension of the model envelope beyond 



this radius, to lower densities and lower temperatures, would 
contribute to a cold, tenuous component. The latter explana- 
tion is favored, since adjusting the temperature profile would 
affect many molecular species, whereas adding low-density ma- 
terial would mostly affect CO and its isotopes. However, two 
test models with extended envelopes, following the same den- 
sity slope as in Sect. A2_ and with a temperature extrapolated 
following r -0 - 4 , suggest that the truncation radius alone cannot 
explain the discrepancy. We define one model with additional 
shells extending to 126kAU and a temperature ranging down 
to 17.0 K, and another extending to 221 kAU and a temperature 
of 13.6 K. The spatial distribution of integrated C 17 3-2 inten- 
sity resulting from the first model is steeper, but still too flat to 
match the observed distribution. The J T^&V profile resulting 
from the second extended model is very similar to the first one, 
suggesting that a further extension of the envelope would not 
bring the modeled intensity distribution closer to the observa- 
tions. Furthermore, we point out that by incorporating the ob- 
served blueshifted component of, a.o., C 17 3-2 emission we 
could be contaminating our view of the quiescent molecular en- 
velope. 

An attempt to reduce line opacity by introducing velocity 
structure (Sect. |4.4| ) is moderately successful in the sense that 
line optical depths are indeed decreased, but not sufficiently to 
explain the spatial distribution of integrated line intensity for 
species such as HCO + and HCN. Moreover, the fact that the 
spectral linewidths produced by the infall model are incompat- 
ible with observations of both optically thick as well as opti- 
cally thin lines (Fig. [9]), could mean that, if any velocity gradient 
exists in the molecular envelope, complete free-fall is not the 
correct approach. Part of the envelope may be radiatively, mag- 
netically or turbulently supported against collapse, while other 
(outer) parts are in fact still collapsing. 

Finally, the flatt ening of the spherical cloud in one spatial 
direction (Sect. 4.5 ) also reduces optical depths, but only when 
viewed under a moderate inclination angle (-20°, see Fig. [12]). A 
high inclination angle has the additional effect of producing spa- 
tial profiles which are more asymmetric than what is observed. 
Hence, if the envelope of AFGL2591 is flattened in reality, we 
expect the system to be viewed almost f ace-o n. We note that both 
trial HCO + abundance values in Sect. 14.51 are consistent with 
lxlO" 8 to within the unc ertainty margin of a factor 2 quoted by 
|Van der Tak etaLl ( [T999l ). 

In addition, it is apparent from Fig. [9] that neither the infall 
model nor the flattened model result in velocity profiles of the 
emission lines which are as narrow as the observed linewidths. 
We note that the discrepancy of the velocity profile for both the 
optically thin and thick lines is larger for the infall model than 
for the static flattened model. It is an interesting result that the 
change from static spherical to static ellipsoidal leaves the line 
shape of the optically thin H 13 CO + 4-3 intact, while it also man- 
ages to let more of the optically thick H 12 CO + 4-3 line intensity 
escape the cloud. 



For HCN specifically, Boonm an et al. ( |2001| ) and Stauber 
|et al.| ( 2005] ) have suggested an abundance which is enhanced 
by as much as a factor 100 in the central few hundred AU. We 
introduce an abundance jump from lxlO -8 at r gas < 230 K to 
lxlO -6 above 230 K in the framework of our static spherical 
model (Sect. |4.2| ), which originally has a constant abundance. 
The thresh old for the abundance jump is motivated by ga s-phase 
chemistry ( Boon man et al.|2 001 ; Va n der Tak et al.|1999| ), rather 
than by sublimation from grain surfaces, which puts the evap- 
oration threshold at -100 K for H 2 ices and at -25 K for CO 
ices. While this abundance jump is in principle a candidate to 
explain the missing line intensity from the central region of the 
envelope (cf. Fig. [TO]), the abundance jump model underproduces 
the observed HCN 4-3 intensity toward the central position by 
~40%, while simultaneously overproducing the off-center inten- 
sity. This is caused by the fact that only the central few hun- 
dred AU enjoys a gas temperature above 230 K, the effects of 
which are washed out by the -15" (15 000 AU) FWHM of the 
telescope beam. Hence, the very inner part of the envelope af- 
fected by the enhanced abundance is too small to modify the 
J r m bdV profile on spatial scales probed by the SLS observa- 
tions. Alternatively, the lack of modeled HCN 4-3 from the cen- 
ter of the envelope can be explained (partly) by the critical den- 
sity, which is on the order of 10 8 cm -3 for HCN 4-3 (at 100 K) 
and only ~10 7 cm -3 for HCO + 4-3. Compared to the maximum 
density of 1.2xl0 7 cm -3 in the central shell of the spherical mod- 
els (Sects. |4~2] |44| ), it is found that HCN is not thermalized. 

The dust continuum modeling in Sect. |4.1| indicates that a 
spherical description of the envelope with a single power law 
density profile cannot simultaneously explain the 450 jim and the 
850 /mi continuum maps. This could indicate a deviation from 
spherical symmetry. Test runs of line radiative transfer models 
with density indices steeper than 1.0 (Sect. |4.3| ) demonstrate that 
model descriptions with a = 1.5 or 2.0 do not provide a more 
satisfying match to the observed emission in spatial and spectral 
dimensions simultaneously. Other envelopes around young high- 
mass stars are generally well fit by power law density profiles 
with a between 1.0 and 2.0, as studied in the samples by | Van der] 
Tak et al.| ( [2000] ) and |Hatchell & Van der Tak] p003] ). Density 
slopes in low-m ass p rotostellar envelopes are found to occupy 
the same range (J0rgens en et al.|[2002 ). The slope of a = 1.0 
for AFGL2591 is on the shallow end of this interval, but a di- 
rect correlation between a an d evolutionary stage is n ot found in 
comparitive studies by, e.g., | Van der Tak et~aL] ([2000). 

Concluding, the model descriptions employed in this study 
only match the observed molecular line intensity distributions 
if the line under consideration is optically thin. From a model- 
ing point of view, signatures from these lines are not sensitive 
to the detailed distribution of the material in the envelope, as 
long as roughly the right amount of material is present at the 
right density and temperature. Conversely, we find a profound 
mismatch for optically thick lines, exactly those that are obser- 
vational probes of the geometry of the envelope and the distribu- 
tion of the material. Below we propose several solutions to the 
shortcomings of our models. 

- Instead of choosing between either velocity structure or non- 
spherical morphology, a combination of these two ingredi- 
ents might yield a closer match of the observed molecular 
emission in the spatial domain. In the spectral domain, how- 
ever, we do not expect a non- spherical inf ailing model to ex- 
hibit line shapes that will match the observations (cf. Fig. [9]). 

- A second route to lower optical depths is the adoption of an 
outflow cavity in the envelope, in combination with a favor- 
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able viewing angle, as p roposed by Van der Ta k et al. (1999]) 
and later investigated by Preibisch et al. '( 2003| ) and Bruderer 
EtaEld2010al. 



10"=10 4 AU 
N 



- A third option is to introduce inhomogeneous ('clumpy') 
structure to alleviate optical depth effects along selected lines 
of sight ( |Wang et al.|[T993l |Spaans & Van Dishoeck|[T997l ). 
The spatial scales of these inhomogeneities should then be 
such that they are largely unresolved by the observations pre- 
sented here, i.e., smaller than ~10 4 AU. 

5.2. Further evidence of substructure 

The models described in Sect. [4] are an attempt to explain 
the large-scale morphology of the molecular envelope of 
AFGL2591 by focusing on a selection of molecular species, but 
we have also not ed ad ditional irregular structure in other tracers 
in Sects. 3.2 and 3.3 Various maps of molecular species which 
we did not model explicitly are therefore interesting to approach 
from a phenomenological perspective. 

The northeastern 'arm' together with the southern indenta- 
tion seen in six different species might indicate a quadrupolar 
outflow. This could be interpreted as a set of two bipolar out- 
flows originating from two separate protostars, instead of the 
previously assumed bipolar outflow associated with a single pro- 
tostar. However, apart from the velocity separation seen in N2H + , 
there is no overwhelming kinematic evidence for this additional 
outflow direction. Alternatively, the 'arms' at position angles 
-130° and -310° (N2H + velocity map in Fig. [6]) may also be 
inf ailing gas streams, perhaps belonging to a larger scale fila- 
mentary structure (Kle ssen et al.|2005"] ). 

The northeastern plume is especially bright in several 
CH 3 OH transitions with varying E up . This is an indication that 
its origin does not lie in excitation effects. Combined with the 
observation that the plume is absent in C2H and the sulfur- 
bearing species, we suggest that chemical effects are at play. 
Grain surface chemistry is known to be an important formation 
route for CH 3 OH ( [Charnley et al.|1997|[Herbst & Van Dishoeck 



2009 ), which could have evaporated from grain mantles under 
the influence of a local shock front. 

The shape and direction (position angle -230°) of 
blueshifted components of 13 CO, C 17 0, HCO + , o-H 2 CO, and 
N 2 H + coincide with the direction of the known approaching 
molecular outflow (Hasegawa & MitchelT| [l995| ). It is interest- 
ing to note that, contrary to most molecules in Figs. [6] and [7] 
where both the envelope and the outflow contribute to the ob- 
served emission, N 2 H + 4-3 in particular appears to trace the 
outflow rather than the bulk envelope. We base this on the sig- 
nificantly flattened morphology of its integrated line emission 
map as well as the velocity gradient which is more pronounced 
for N 2 H + than for the other molecules observed here. The N 2 H + 
molecule has previously been used as a tracer of dense mate- 
rial (e.g.,|Caselli et al.|2002a|b[|Crapsi et al.|2004[|Di Francesco 



et al.|2004 ), but usually from observations of the 7=1-0 or 3-2 



transition. We not e that maps of N 2 H + 4-3 are rare in the lit- 
erature. Moreover, Busq uet et al.| ( |2011| ) have found significant 
depletion of N 2 H + w.r.t. NH 3 in dense parts of hot cores, again 
derived from N 2 H + 1-0 observations. So apart from the pro- 
nounced velocity structure, density dependent depletion could 
contribute to the atypical distribution of N 2 H + emission seen in 
our map (Fig. [2]). Conve rsely, the N 2 H + emission could trace 
the evaporation of N 2 (Bergin et al. 2002 ) along outflow cavity 
walls. 

Although only marginally significant w.r.t. the pointing un- 
certainties, the fact that high-density tracers such as HNC and 




warm methanol 
spots 



Fig. 14. Overview of observed substructure in the circumstellar 
envelope of AFGL2591. The 'outer' and 'inner' envelopes, as 
well as the four- arm structure of N 2 H + , the warm CH3OH spots, 
and the CH3OH plume are described in this work. The directions 
of the approaching (blue) and receding (red) outflows are indi- 
cated. The yellow star marks the position of the central heating 
source (infrared source from |Van der Tak et al.|1999]). VLA1 and 
VLA2 are radio continuum sources detected by Trinidad et al. 
(2003]). The near-IR loops from Preibis ch" et al.f ( |2003| ) are h> 
dicated in dashed yellow. For reference, the fields of view of 
JCMT/HARP-B and the interferometric maps from VLA and 
OVRO are indicated. 



CN have peak positions offset to the south w.r.t. other species 
could be a result of density anisotropies. C 2 H being offset to 
the northwest could arise from chemical processes influenced 
by non-isotropic UV-radiation. 

The warm CH 3 OH 12i-12 A T transition, with £ up =197K, 
shows a two-peak structure (see Sect. [3]), indicative of two sepa- 
rated heating sources. While four emission sources with separa- 
tions < 6" have been seen in interferometric observations of this 
sourc e ( [Van der Tak et al.|1999HTriri idad et al. 2003 ; Benz et al.| 
|2007| ), the warm components separated by > 10" (>10000 AU) 
described here have not been seen before. Although the 12i-12o 
transition is not listed in studies of masing CH3OH lines (e.g., 
|Cragg et al.|1 992), we cannot exclude the possibility that it is a 
maser. If it is, then a line width narrower than the 5 km s" 1 that we 
observe would need to be confirmed by higher angular resolution 
observations. An intrinsically narrow line profile is potentially 
smeared out by our 15" beam, which encompasses emission 
from a considerable portion of the molecular cloud. Moreover, 
excitation analysis of an ensemble of CH 3 OH lines, planned in 
Van der Wiel et al. (2011b, in prep.) will reveal whether or not 
the 12i— 12q emission fits in the picture of thermal emission. 



5.3. Conclusions and outlook 

To summarize, this study uses the unbiased spectral imaging 
from the JCMT Spectral Legacy Survey to investigate the struc- 
ture of the protostellar envelope of AFGL2591. We employ 35 
spatially resolved molecular line signatures to investigate the en- 
velope in three dimensions: two spatial axes with a resolution 
element of ~ 10 4 AU in projected distance, and the line-of-sight 
velocity axis with a resolution of 1 km s" 1 . 
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Figure 14 provides an overview of selected observations 
from the literature and from this study which are indicative of 
substructure in the envelope of AFGL2591. With this overview 
we aim to support the suggestion that inhomogeneous substruc- 
ture at scales of a few 1000AU and smaller should be present 
in molecular envelopes around protostars. The two separated 
spots of warm CH3OH, a warped asymmetric plume feature, 
a flattened morphology of warm SO2, and the multiple 'arm' 



structures seen in CO, HCN, H 2 CO, N 2 H + , and CS (Sect. [5^2 ) 
could arise from local underdense patches with associated lower 
cumulative column densities, higher (UV)-flux conditions, and 
consequently higher gas temperatures and altered chemistry. At 
the same time, such a picture could resolve optical depth ef fects 
for the quiescent envelope tracers discussed in Sects. [4] and 5.1 



Therefore, our results warrant further theoretical modeling ef- 
forts for protostellar envelopes in the direction of clumpy models 
and of flattened model structures, possibly including cavities and 
velocity structure. We emphasize the value of the spectral imag- 
ing survey provided by the SLS, which allows for the simultane- 
ous employment of several molecular tracers to characterize the 
physical properties of the large-scale envelope of AFGL2591. 

An unbiased Herschel/HIFI spectral survey of AFGL2591 in 
the 490-1240 GHz regime is currently under analysis (Van der 
Wiel et al. 2011a, in prep.). This will give insights into the ex- 
citation conditions in the molecular envelope, but does not con- 
tain any spatial information. Moreover, spectral mapping of a 
consistent set of molecular tracers such as those presented here, 
but at higher angular resolution, will be needed to confirm the 
'clumpy' picture sketched above. The set of tracers should be 
chosen to be sensitive to the overall quiescent envelope as well 
as to pockets of warm (> 200 K) actively heated gas. Although 
our specific target is well into the northern sky, the obvious can- 
didate observatory to study similar high-mass star-forming en- 
velopes in the southern sky is the Atacama Large Millimeter and 
Submillimeter Arra}^] This facility will be able to map high- 
/ lines of simple molecules like CO, HCN, HCO + - and their 
rarer isotopic variants - at subarcsecond angular resolution. 
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